Regressão Linear Múltipla em R, parte II
Reinaldo Zezela, DMI, Universidade Eduardo Mondlane, Mozambique
11 March 2020
Photo by Isaac Smith on Unsplash
1.0 Regressão linear Múltipla
A regressão linear simples é uma abordagem útil para prever uma resposta com base em uma única variável preditora. No entanto, na prática, muitas vezes temos mais de um preditor. Por exemplo, nos dados de publicidade visto no capítulo anterior (Regressão Linear Simples), examinamos a relação entre vendas e o gasto em publicidade na TV. O dados também mostram o montante em dinheiro gasto em publicidade na rádio e jornal, então estamos interessados em saber se as vendas estão associadas a uma dessas duas mídias.
Portanto, como podemos estender a nossa análise dos dados de publicidade para acomodar esses dois preditores adicionais?
Em geral, suponha que temos \(p\) preditores distintos. Então, o modelo de regressão linear múltipla assume a forma \[ Y = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2} + ... + \beta_p X_{p} + \epsilon \], onde \(X_{j}\) representa o \(j\)-ésimo predictor e \(\beta_j\) quantifica a associação entre essa variável e a resposta. Nós interpretamos \(\beta_j\) como o efeito médio em \(Y\) de um aumento de uma unidade em \(X_{j}\), mantendo todos os outros preditores fixos. No presente caso, sobre os dados de publicidade temos a seguinte forma: \[ Y = \beta_0 + \beta_1*TV + \beta_2* radio + \beta_3 * newspaper + \epsilon \].
1.1 Estimar os coeficientes de regressão
Como foi o caso no ajuste de regressão linear simples, os coeficientes de regressão \(β0, β1,. . . , βp\) são desconhecidos e devem ser estimados. Dadas as estimativas \(\hat{\beta_0}, \hat{\beta_1},..., \hat{\beta_p}\), podemos fazer previsões usando a fórmula \[ \hat{y} = \hat{\beta_0} + \hat{\beta_1} X_{1} + \hat{\beta_2} X_{2} + ... + \hat{\beta_p} X_{p} + \epsilon \].
Os parâmetros são estimados usando a mesma abordagem de mínimos quadrados que vimos no contexto da regressão linear simples. Nós escolhemos \(β0, β1,. . . , βp\) para minimizar a soma dos quadrados dos resíduos \[ RSS = \sum_{i=1}^{n}{(y_i - \hat{y_i})^2} \].
Em código R fazemos de seguinte modo:
1o) carregar os dados à semelhança que fizemos na Regressão Linear Simples;
2o) usar a fórmula lm(). Podemos executar ?lm para ajuda.
# importar os pacotes requeridos
library(tidyverse) # para manipulação de dados## Warning: package 'tidyverse' was built under R version 3.6.3## -- Attaching packages ------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.0 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.4
## v tidyr 1.0.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## Warning: package 'tibble' was built under R version 3.6.2## Warning: package 'tidyr' was built under R version 3.6.3## Warning: package 'purrr' was built under R version 3.6.3## Warning: package 'dplyr' was built under R version 3.6.3## -- Conflicts ---------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(pander) # para formatação em R Markdown## Warning: package 'pander' was built under R version 3.6.2# carregar os dados
publicidade_dados <- read.csv(file = 'data/Advertising.csv')Como sempre temos que examinar os dados.
publicidade_dados %>% glimpse()## Observations: 200
## Variables: 5
## $ X <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
## $ TV <dbl> 230.1, 44.5, 17.2, 151.5, 180.8, 8.7, 57.5, 120.2, 8.6, 1...
## $ radio <dbl> 37.8, 39.3, 45.9, 41.3, 10.8, 48.9, 32.8, 19.6, 2.1, 2.6,...
## $ newspaper <dbl> 69.2, 45.1, 69.3, 58.5, 58.4, 75.0, 23.5, 11.6, 1.0, 21.2...
## $ sales <dbl> 22.1, 10.4, 9.3, 18.5, 12.9, 7.2, 11.8, 13.2, 4.8, 10.6, ...O output mostra variável \(X\) que não é relevante para o nosso estudo, portanto vamos remove-lá.
publicidade_tbl <- publicidade_dados[ , -1]
publicidade_tbl %>% head() %>% pander()| TV | radio | newspaper | sales |
|---|---|---|---|
| 230.1 | 37.8 | 69.2 | 22.1 |
| 44.5 | 39.3 | 45.1 | 10.4 |
| 17.2 | 45.9 | 69.3 | 9.3 |
| 151.5 | 41.3 | 58.5 | 18.5 |
| 180.8 | 10.8 | 58.4 | 12.9 |
| 8.7 | 48.9 | 75 | 7.2 |
Para determinar os coeficientes de regressão recomenda-se:
1o passo: sempre visualize os seus dados.
plot(publicidade_tbl)
O output mostra o seguinte:
- a representação gráfica entre as variáveis. Podemos notar o mesmo gráfico entre vendas \(sales\) e gastos em \(TV\) discutido na aula anterior, Regressão Linear Simples.
- o gráfico entre vendas e rádio, portanto \(sales\) e \(radio\) mostra um relacionamento aparentemente linear, enquanto que entre o jornal \(newspaper\) e as vendas \(sales\) mostra claramente um relacionamento não linear.
2o passo: usar a fórmula.
multiple_model <- lm(sales ~ TV + radio + newspaper, data = publicidade_tbl)
summary(multiple_model)##
## Call:
## lm(formula = sales ~ TV + radio + newspaper, data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8277 -0.8908 0.2418 1.1893 2.8292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.938889 0.311908 9.422 <2e-16 ***
## TV 0.045765 0.001395 32.809 <2e-16 ***
## radio 0.188530 0.008611 21.893 <2e-16 ***
## newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-161.2 Interpretação
A primeira etapa na interpretação da análise de regressão múltipla é examinar F-statistic e o p-value associado, na parte inferior do output do resumo do modelo.
Em nosso exemplo, pode-se ver que o p-value de F-statistic é \(<2,2e-16\), o que é altamente significativo. Isso significa que, pelo menos, uma das variáveis preditivas está significativamente relacionada à variável de resultado.
Para averiguar quais variáveis preditivas são significativas, é possível examinar a tabela de coeficientes, que mostra a estimativa dos coeficientes beta de regressão e os valores de t-statitic e p-values associados:
summary(multiple_model)$coefficients %>% pander()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.939 | 0.3119 | 9.422 | 1.267e-17 |
| TV | 0.04576 | 0.001395 | 32.81 | 1.51e-81 |
| radio | 0.1885 | 0.008611 | 21.89 | 1.505e-54 |
| newspaper | -0.001037 | 0.005871 | -0.1767 | 0.8599 |
O output mostra o seguinte:
- Para um dado preditor, o p-values avalia se existe ou não associação significativa entre o preditor e a variável de resultado, ou seja, se o coeficiente beta do preditor é significativamente diferente de zero.
Pode-se observar que a alteração no orçamento de publicidade de \(TV\) e da \(radio\) está significativamente associada a alterações nas vendas \(sales\), enquanto as alterações no orçamento do jornal \(newspaper\) não estão significativamente associadas às vendas \(sales\), visto que o seu p-values \((0.8599)\) é maior que \(0.05\).
Para uma dada variável preditora, o coeficiente beta pode ser interpretado como o efeito médio em \(y\) de um aumento de uma unidade no preditivo, mantendo todos os outros preditores fixos.
Por exemplo, mantendo os gastos fixos em publicidade na \(TV\) e jornal \(newspaper\), e gastar um valor adicional de \(1,000.00\) US dolar em publicidade na \(radio\) leva a um aumento nas vendas em aproximadamente \(0,1885 * 1000 = 189\) unidades de venda, em média.
O coeficiente de \(TV\) sugere que, para cada aumento de \(1,000.00\) dolares no orçamento de publicidade em \(TV\), mantendo todos os outros preditores constantes, podemos esperar um aumento de \(0,045 * 1000 = 45\) unidades de vendas, em média.
Descobrimos que o jornal \(newspaper\) não é significativo no modelo de regressão múltipla. Isso significa que, para uma quantia fixa de orçamento de publicidade em \(TV\) e em jornal \(newspaper\), as alterações no orçamento de publicidade de jornal \(newspaper\) não afectarão significativamente as unidades de vendas.
Portanto, como a variável jornal \(newspaper\) não é significativa, é possível removê-la do modelo:
multiple_model_v2 <- lm(sales ~ TV + radio, data = publicidade_tbl)
summary(multiple_model_v2)##
## Call:
## lm(formula = sales ~ TV + radio, data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7977 -0.8752 0.2422 1.1708 2.8328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.92110 0.29449 9.919 <2e-16 ***
## TV 0.04575 0.00139 32.909 <2e-16 ***
## radio 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Finalmente, nossa equação de modelo pode ser escrita da seguinte forma: \[sales = 2.921 + 0.046 * TV + 0.188 * radio\].
O intervalo de confiança do coeficiente do modelo pode ser extraído da seguinte forma:
confint(multiple_model, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 2.32376228 3.55401646
## TV 0.04301371 0.04851558
## radio 0.17154745 0.20551259
## newspaper -0.01261595 0.01054097Também podemos obter os intervalos de confiança em torno dessas estimativas de coeficiente, como fizemos no modelo de regressão linear simples. Aqui vemos como o intervalo de confiança para o jornal \(newspaper\) inclui o valor zero \(0\), o que sugere que não podemos assumir que a estimativa de coeficiente de \(-0.001037\) seja diferente de \(0\).
1.3 Avaliação da precisão do modelo
Como vimos na regressão linear simples, a qualidade geral do modelo pode ser avaliada examinando o coeficiente de determinação R-squared (\(R^2\)) e o Erro Residual Padrão (RSE).
Coeficiente de Determinação: \(R^2\)
Na regressão linear múltipla, o \(R^2\) representa o coeficiente de correlação entre os valores observados da variável de resultado (\(y\)) e os valores ajustados (isto é, previstos) de \(y\). Por esse motivo, o valor de \(R^2\) sempre será positivo e variará de zero a um.
\(R^2\) representa a proporção de variação, na variável de resultado \(y\), que pode ser prevista pelo conhecimento do valor das variáveis \(x\). Em outras palavras, é a proporção de variação de \(y\) que pode ser explicada pela equação de regressão. Um valor de \(R^2\) próximo a 1 indica que o modelo explica grande parte da variação na variável de resultado.
Um problema com o \(R^2\) é que ele sempre aumentará quanto mais variáveis forem adicionadas ao modelo, mesmo que essas variáveis estejam fracamente associadas à resposta (James et al. 2014). Uma solução é ajustar o \(R^2\) levando em consideração o número de variáveis preditoras.
O ajuste no valor de \(R^2\) “Adjusted R-squared”, conforme mostra no output do resumo é uma correção para o número de variáveis \(x\) incluídas no modelo de previsão.
No nosso exemplo, com variáveis preditivas de \(TV\) e \(radio\), o \(R^2\) ajustado (Adjusted R-squared: 0.8962), significando que \(89.62\)% da variação na medida de vendas pode ser prevista pelos orçamentos de publicidade em \(TV\) e \(radio\).
Este modelo é melhor que o modelo linear simples, apenas com \(TV\) discutido na aula anterior (regressão linear simples), que tinha um \(R^2\) ajustado de \(61\)%.**
Erro padrão residual (RSE) ou sigma:
A estimativa do RSE fornece uma medida de erro de previsão. Quanto menor for o RSE, mais preciso é o modelo (nos dados em mãos).
A taxa de erro pode ser estimada dividindo o RSE pela variável de resultado médio:
round(sigma(multiple_model_v2)/mean(publicidade_tbl$sales), 3)## [1] 0.12No nosso exemplo de regressão múltipla, o RSE é \(1.681\), correspondendo a uma taxa de erro de \(12\)%.
Novamente, isso é melhor que o modelo simples, com apenas a variável \(TV\), onde o RSE era de \(3.259\) (~ \(23\)% de taxa de erro) discutido na aula anterior (regressão linear simples).
2.0 Diagnóstico
Além de observar as estatísticas do \(RSE\) e \(R^2\) que acabamos de discutir, pode ser útil visualizar os dados. Resumos gráficos podem revelar problemas com um modelo que não são visíveis a partir de estatísticas numéricas.
library(ggfortify)## Warning: package 'ggfortify' was built under R version 3.6.3autoplot(multiple_model_v2)
Observe que há um padrão claro nos resíduos, o que põe em causa o pressuposto sobre a linearidade. O gráfico sobre a normalidade também é preocupante. Portanto, esse padrão não linear pronunciado não pode ser modelado com precisão usando a regressão linear. Sugere um efeito de sinergia ou interação entre os mídias, em que combinar a mídia resulta em um maior impulso para as vendas do que usar qualquer meio único.
Na próxima secção, discutiremos a extensão do modelo linear para acomodar tais efeitos sinérgicos através do uso de termos de interação. Outro aspecto importante para lidar com um modelo não linear será proceder a transformação das variáveis para função logarítimica. Matéria que será objecto de análise nos próximos capítulos.
3.0 Interação
A equação de regressão linear múltipla, com efeitos de interação entre dois preditores (\(x_1\) e \(x_2\)), pode ser escrita da seguinte forma: \[ y = b_0 + b_1 * x_1 + b_2 * x_2 + b_3 * (x_1 * x_2)\]
Considerando o nosso exemplo, ele se torna:\[ sales = b_0 + b_1 * TV + b_2 * radio + b_3 * (TV * radio)\]
\(b_3\) pode ser interpretado como o aumento da eficácia da publicidade na \(TV\) para um aumento de uma unidade na publicidade no \(radio\) (ou vice-versa).
No R, a interação entre as variáveis é usando o operador ’*’:
# Build the model
# Use this:
interaction_model <- lm(sales ~ TV + radio + TV:radio, data = publicidade_tbl)
# Or simply, use this:
interaction_model <- lm(sales ~ TV*radio, data = publicidade_tbl)
# Summarize the model
summary(interaction_model)##
## Call:
## lm(formula = sales ~ TV * radio, data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3366 -0.4028 0.1831 0.5948 1.5246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.750e+00 2.479e-01 27.233 <2e-16 ***
## TV 1.910e-02 1.504e-03 12.699 <2e-16 ***
## radio 2.886e-02 8.905e-03 3.241 0.0014 **
## TV:radio 1.086e-03 5.242e-05 20.727 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared: 0.9678, Adjusted R-squared: 0.9673
## F-statistic: 1963 on 3 and 196 DF, p-value: < 2.2e-16Interpretação:
- Pode-se observar que todos os coeficientes, incluindo o coeficiente do termo de interação, são estatisticamente significativos, sugerindo que existe uma relação de interação entre as duas variáveis preditivas (publicidade em \(TV\) e \(radio\)).
Nossa equação de modelo fica assim: \[ sales = 6.75 + 0.019 * TV + 0,029 * radio + 0,0009 * TV * radio \]
Podemos interpretar isso como um aumento de \(1,000.00\) dólares em publicidade na \(TV\) associado a um aumento nas vendas de \((b_1 + b_3 * radio) * 1000 = 19 + 0,9 * unidadesdeRadio\). E um aumento na publicidade na \(radio\) de \(1,000.00\) dólares será associado a um aumento nas vendas de \((b2 + b3 * TV) * 1000 = 29 + 0,9 * unidades de TV\).
Nota: Observe que, algumas vezes, o termo interação é significativo, mas não são os principais efeitos. O princípio hierárquico afirma que, se incluirmos uma interação em um modelo, também devemos incluir os efeitos principais, mesmo que os valores de p-value associados a seus coeficientes não sejam significativos (James et al. 2014).
Fazendo uma comparação entre os modelos aditivos e de interação podemos constatar que:
o Erro padrão residual (RSE) reduziu de (Residual standard error: 1.681) - modelo aditivo para (Residual standard error: 0.9435) - modelo interactivo. Ou seja, quanto menor fôr o RSE melhor os dados se ajustam ao modelo.
Além disso, o valor de \(R^2\) do modelo de interação é de \(97\)% em comparação com apenas \(90\)% no modelo de aditivo.
Esses resultados sugerem que o modelo com o termo de interação é melhor que o modelo que contém apenas efeitos principais. Portanto, para esses dados específicos, devemos optar pelo modelo de interação para predição.
4.0 Discussão
Este capítulo descreve o modelo de regressão linear múltipla.
Observe que, se tiver muitas variáveis preditivas em seus dados, não precisará necessariamente digitar o nome deles ao calcular o modelo.
Para calcular a regressão múltipla usando todos os preditores no conjunto de dados, basta digitar o seguinte:
multiple_model_v3 <- lm(sales ~ ., data = publicidade_tbl)
summary(multiple_model_v3)##
## Call:
## lm(formula = sales ~ ., data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8277 -0.8908 0.2418 1.1893 2.8292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.938889 0.311908 9.422 <2e-16 ***
## TV 0.045765 0.001395 32.809 <2e-16 ***
## radio 0.188530 0.008611 21.893 <2e-16 ***
## newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16Se desejar executar a regressão usando todas as variáveis, excepto uma, digamos jornal \(newspaper\), digite o seguinte:
multiple_model_v4 <- lm(sales ~ . -newspaper, data = publicidade_tbl)
summary(multiple_model_v4)##
## Call:
## lm(formula = sales ~ . - newspaper, data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7977 -0.8752 0.2422 1.1708 2.8328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.92110 0.29449 9.919 <2e-16 ***
## TV 0.04575 0.00139 32.909 <2e-16 ***
## radio 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Como alternativa, pode usar a função update():
multiple_model_v5 <- update(multiple_model_v3, ~. -newspaper)
summary(multiple_model_v5)##
## Call:
## lm(formula = sales ~ TV + radio, data = publicidade_tbl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7977 -0.8752 0.2422 1.1708 2.8328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.92110 0.29449 9.919 <2e-16 ***
## TV 0.04575 0.00139 32.909 <2e-16 ***
## radio 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Este capítulo descreve também como calcular a regressão linear múltipla com efeitos de interação. Os termos de interação devem ser incluídos no modelo se forem significativos.
O processo de transformação das variáveis e avaliação de seus impactos serão discutidos no próximo capítulo.
5.0 Referências
James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2014. An Introduction to Statistical Learning: With Applications in R. Springer Publishing Company, Incorporated.
Bruce, Peter, and Andrew Bruce. 2017. Practical Statistics for Data Scientists. O’Reilly Media.